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ABSTRACT 



A molecular dynamics computer simulation was used to investigate several 
techniques of generating liquid Cu targets. The target with the best liquid 
characteristics was subjected to one KeV, argon ion bombardment as a preliminary 
study of the sputtering of liquids. The techniques of warming by impulse and 
warming by initially displacing atoms from their equilibrium positions were 
compared. Both methods produced targets with good liquid properties. The energy 
became equally partitioned between kinetic and potential energy and all targets 
equilibrated within 400 fs. The range of a typical atom during the time of 
equilibration was found to be restricted to its initial neighborhood. The preliminary 
sputtering study resulted in a sputtering yield increase of 40% over the solid target, 
for a low index crystal plane. 
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I. INTRODUCTION 



A. HISTORICAL OVERVIEW OF SPUTTERING 

In 1851, Plucker [ref. 1] observed that the gas in x-rays tubes was continually 
removed. He attributed this phenomenon to the ionization of residual gases and 
their consequent absorption by the inner surface of the x— ray tube. In 1852, Grove 
[ref. 2] noticed that the x— ray tube surface struck by these particle ions was eroded 
due to the removal of target material. He called this phenomena cathode sputtering. 

At that point in history, the action of target erosion or sputtering was 
considered undesirable in lab equipment, and was only interesting so long as one 
could minimize its unwanted effects. Consequently, sputtering was not investigated 
systematically for more than fifty years. 

About fifty years after Grove reported his cathode sputtering findings, 
Goldstein [ref. 3] presented conclusive evidence that sputtering was indeed caused 
by the positive atoms of the discharge impacting on the cathode target. 

In 1908, Stark [ref. 4] advanced the concept of the individual sputtering event 
on an atomic scale. He developed a collision theory which treated sputtering as a 
sequence of binary collisions initiated by the collision of one incident ion at a time. 
He also presented a second theory called the hot-spot model which attributed the 
sputtering action to localized high temperature heating of the target and 
evaporation of atoms. In 1921, Thompson [ref. 5] proposed that sputtering was 
caused by the release of radiation as the bombarding ions struck the target. The 
following year Bush and Smith [ref.6] suggested that sputtering was caused by the 
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expansion of gas adsorbed by the target material, and Kingdon and Langmuir 
[refs. 7,8] conducted a sputtering experiment which yielded ejecta from the surface 
layers of the target. They bombarded thoriated tungsten with ions in a glow 
discharge tube. The experiment demonstrated that the bulk of the ejecta was 
coming from the thin film of thorium on the target's surface instead of the tungsten 
substrate. Their results suggested a momentum transfer ejection mechanism for 
sputtering. 

In 1926, Von Hippie and Blechschmidt [refs. 9—12], proposed a theory that 
described sputtering as an evaporation of the surface atoms. Von Hippie extended 
Stark's hot-spot model and attempted to formulate a sputtering theory on the basis 
of local heating. He expressed the view that local heating was the only feasible way 
to explicitly treat the statistics of the complex collisions occurring in a sputtering 
event. 

Approximately ten years later, Lamar and Compton [ref. 13] published A 
Special Theory of Cathode Sputtering which led to the thermal spike concept. The 
thermal spike was based on momentum transfer between the incomming ion and the 
lattice atoms. This theory suggested that a long-lived high temperature volume 
persisted in the target even after the collision cascade was completed. 

In 1931. Guntherschulze and Meyer [refs. 14,15] were the first to recognized the 
importance of minimizing ambient pressure in sputtering experiments. They used a 
high vacuum tube and took the precaution of removing several top layers from the 
target to clean its surface. Consequently, their results were the first to satisfy the 
conditions for reproducible sputtering yield determination. About ten years later, 



2 



Penning and Mobious [ref. 16] conclusively demonstrated the significant effects of 
ambient pressure on the sputtering yield. 

By the early 1950’s, a renewed interest in sputtering coalesced in the scientific 
community. This was brought about primarily by the fact that sputtering was 
found to have technological value. It was discovered that the ion— trapping process 
in sputtering could be used as a pumping effect for low pressure electronic systems. 
It was also found that sputtering could be used to clean target surfaces. Sputtering 
became useful to industry and its status was raised from that of a laboratory 
hindrance to one meriting serious scientific research effort. 

In 1952, Key well [ref. 17] formulated a sputtering theory which made use of 
existing neutron transport theories originally developed for nuclear reactor design. 
Keywell's work, as well as subsequent calculations by Harrison [ref. 18] made 
important contributions to the field of sputtering by introducing probability 
concepts in terms of collisional cross sections. 

In 1954, Wehner [ref. 19] published his findings on crystal structure effects in 
the flux of sputtered particles. This was a significant advance in that it proved 
conclusively that Stark's hot-spot model was incapable of fully explaining the 
sputtering mechanism. Wehner's conclusions fortified collisional theory as an 
important part of sputtering theory. His findings created great interest in the 
collisional aspect of sputtering and served to point the way for following research 
efforts. 

In the decade that followed, (mid 50's to mid 60's) the main emphasis in 
sputtering research w r as directed towards the study of crystal lattice effects. Such 
studies produced two main theories. These were the channeling theory and the 
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focusing collision theory. The channeling theory was successful in describing the 
angular variation of the yield when the incident ions where closely aligned to one of 
the target's principal axis or planes. The theory failed for general alignment of the 
incident ions. The focusing collision theory was the product of Silsbee [ref. 20]. His 
theory showed that momentum was transferred within the target's crystal structure 
along preferred directions. He introduced the idea of momentum transport without 
requiring mass displacement. This concept was later refered to as focusons. 
Focusing collision theory was relatively successful in modelling the observed ejection 
patterns from materials with a high degree of symmetry, such as cubic crystals, but 
only if bombarded with high energies. This theory, however, failed to match the 
Wehner spots and proved to be particularly inadequate for targets with low 
symmetry such as hexagonal crystal structures. Experimental results demonstrated 
that the sputtering yield was significantly dependent on the crystallographic 
orientation of the target with respect to the incident ion beam. 

The theory of collisional cascades was further advanced by the work of Sigmund 
[refs. 21,22]. He used Lindhard's range theory [ref. 23] and proposed that for 
amorphous solids, the collision cascades could be described by a Boltzmann 
transport equation [ref. 24]. His equations required the target surface to have a 
disordered structure, not the long straight rows of atoms intersecting the surface 
used in the focusing model. He obtained first order asymptotic solutions for cascade 
events produced by ions of medium to large atomic number with energy in the KeV. 
His theory became the reference standard for sputtering yield measurements despite 
its limitation of application to polycrystalline solids with randomly oriented 
crystallites as an amorphous approximation. This idea was simultaneously 
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investigated by Thompson [ref. 25]. He proposed that the ejected atom was affected 
by the surface attraction of the target's surface, causing a deviation on the velocity 
vector of the outgoing particle. This resulted in a distortion of the angular 
distribution of ejected particles. Another important idea developed from collisional 
cascades was that of replacement collision sequences. In replacement collision 
sequences, the moving atom replaces an atom on its lattice site. The vacated atom 
the proceeds to strike and replace the next atom in the row. The sequence 
propagates as each atom replaces the next. 

Despite its long history, sputtering researchers have not been able to formulate 
an analytical theory which can fully explain sputtering. There are still differences 
between the theoretical predictions and experimental data, and it seems that an 
analytically simple form, in all likelyhood, will not be sufficiently flexible to 
correctly describe all aspects of sputtering. The crux of the problem lies in the fact 
that in many cases, the theory requires the solution of a many— body problem with 
multiple interactions. This is a formidable task and one which historically has 
yielded, at best, only approximate solutions. Another approach was required which 
could explain and predict sputtering events based on a few simple laws. With the 
advent of the high speed coputer capable of handling the voluminous amount of 
required calculations, some researchers turned to computer simulations. 

B. COMPUTER SIMULATIONS OF SPUTTERING 

1. Historical Overview 

By the late 1960's, the once curious laboratory observation of Plucker 
[ref. 1] and Grove [ref. 2], had been studied for approximately 100 years and 
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sputtering theories had evolved into dependable sputtering yield predictors. 
Especially successful in determining sputtering yield were the statistical theories of 
Thompson [ref. 25], ejected atom energy distribution function dN/dE, and Sigmund 
[ref. 24], sputtering yield as a function of energy, Y(E). 

The advent of the high speed computer made it practical to perform the 
many calculations required to treat the collision cascade problem in a three 
dimensional crystal. And so, computer simulation made its debut in sputtering 
research. 

A very important early computer simulation was the work of Gibson, et al 
[ref. 26] in 1960. They developed a computer model which simulated the motion of 
primary knock—on atoms in a copper monocrystal target. Their model assigned an 
arbitrary kinetic energy and direction to the struck atom as a means to simulate the 
collision effect. Atomic interactions were treated as binary collisions and the 
resulting motion obtained through Newtonian mechanics. Their results 
demonstrated the ability of computer simulations to isolate elementary sputtering 
prcesses for investigation. Two years later Robinson and Oen [ref. 27], utilizing a 
similar code, discovered the channeling effect. What makes this discovery so 
important is the fact that channeling had not yet been observed in the laboratory. 
For the first time, a computer simulation had predicted an important dynamic 
sputtering mechanism. Channeling was later verified experimentally. 

In 1967, Harrison, Levy, and Effron [ref. 28] simulated the bombardment 
of a copper monocrystal by argon ions using repulsive pair potential interactions. 
Their findings demonstrated that the bulk of the sputtering yield came from the 
first three layers of the top surface. This model was later refined to include 
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attractive potential functions for the target [refs. 29—30]. This potential which 
includes attractive and repulsive interactions produced a dynamically stable crystal. 
In 1978 Garrison, Wi nograd and Harrison [ref. 31] published the result of a 
comprehensive simulation study of atomic and molecular ejection from a copper 
crystal with adsorbed oxygen atoms. The simulation provided a detailed description 
of the ejection mechanism. The results also allowed them to determine whether 
molecules were ejected as clusters or whether they combined together after leaving 
the target's surface. They also determined the effects of adatom placement on 
molecule formation. 

In its short history computer simulation has successfully predicted and 
explained various specific sputtering events. Results of simulations have been both 
praised and strongly criticized by some theoreticians who argue that simulations 
neglect important factors. Nevertheless, the value of computer simulation has been 
recognized by the overall scientific community and it is now considered an 
important research area in the field of sputtering. 

2. General Concepts and Development 

Most real systems are so complicated that a complete analytical 
description is virtually impossible to achieve, unfortunately, sputtering falls into 
this category of complexity. 

The modelling of a real collision cascade such as a sputtering event is an 
attempt to describe a complex system in terms of idealized and highly simplified 
descriptors. Invariably, many factors must be neglected in order to simplify the 
model the selected set of descriptors must adequately determine the system's 
relevant behavior. Once this is achieved, the result is a simulation of a real system. 



Computer simulations of atomic collisions allow researchers to concentrate 
on the basic physics of a system without the analytic constraints of statistical 
theories. It is a tool that can be used to test the applicability of a theory. 

There are two major approaches used in computer simulations. These are 
the time-step and discrete event models. The discrete event model proceeds from 
event to event skipping over their time separation. It maintains a list of possible 
future events from which it determines the next event. This list is constantly being 
updated. This model works well when the events are sufficiently separated in time. 
In contrast, the time— step model advances the clock for a short interval first, and 
then computes all the interactions that occur in that time, updating all processes 
before the next time step. This model is the obvious choice for systems with 
simultaneous interactions. The time— step logic's chief drawback is that it requires 
significantly more computer running time than the discrete event program logic. 

In sputtering applications, these two approaches lead to the two principal 
methods of simulation; the binary collision (BC) and the multiple interaction (MI) 
program logics. 

The BC method [ref. 32] uses the discrete event model. The basic 
assumption of this type of simulation is that each particle interacts with only one 
other particle (usually assumed stationary) at a time. These simulations are 
inherentelly restricted to be linear calculations. 

The multiple interaction simulation uses the time— step model and 
Newton's equations of motion are numerically solved for many particles. This is the 
method used in this thesis investigation. 
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C. SPUTTERING OF LIQUIDS 

The sputtering of liquid targets has gained a great deal of attention in the 
scientific community. Experimental results of sputtering have shown that solid 
targets may undergo a phase transition to liquid under ion bombardment. 
Therefore, there is a practical as well as scientific interest to understand and be 
able to model the characteristics of a liquid surface. 

1. Summary of Experimental Results 

In 1962, Whener et al [ref. 33] reported their measurements on tin 
sputtered by argon ions. Their results showed a 40% larger yield for liquid tin at an 
ion energy of 0.2 KeV and a 6% smaller yield for an ion energy of 0.4 KeV. This 
was the earliest evidence of an energy dependent yield profile with a maximum 
sputtering yield peak. Eight years later Krutenat and Panzera [ref. 34] conducted a 
similar sputtering experiment on solid and liquid tin. Their results were in general 
agreement with those of Whener, showing that the solid yields below 0.2 KeV is 
twice that of the liquid and that a crossover occurs at 0.375 KeV above which the 
liquid yield is only 15% higher. Their results also implicated surface dynamics in 
sputtering mechanisms. 

In 19S2, Cooper and Hurst [ref. 35] bombarded liquid and solid indium 
with argon ions and discovered that above 0.107 KeV the liquid yield was 10% 
higher than the solid. They noted that the shape of the yield versus ion energy for 
both the liquid and solid were similar. They attributed the higher yield of the 
liquid to its lower binding energy. 

Dumke et al [ref. 36] published their results on the sputtering of a 
gallium— indium eutectic alloy (surface monolayer greater than 94% indium) in the 
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liquid phase. Their findings indicated a ratio indium to gallium sputtering yield 
which was 2S times greater than expected from the stoichiometry. They were able 
to conclude that 85% of the sputtered atoms came from the surface monolayer. 

2. Historical Summary of Liquid Models 

In order to study the effects of sputtering on a liquid surface using 
computer simulations, one first needs to develop a reasonable model of a liquid. 

One very early model of a simple liquid was developed by Bernal in 1960 
[ref. 37], He stipulated that the structure of a simple liquid can be determined 
simply by volume exclusion. He then proposed a zeroth order model in which the 
atoms were considered hard spheres and their local structure determined by the 
restriction that no two atoms can approach each other by less than one atom 
diameter. Early Bernal models were constructed [refs. 38,39] by pouring thousands 
of steel balls into a deformable container (originally a football bladder). The 
ensemble was then bound with rubber strips and kneaded to maximize its density. 
The contents were then fixed by pouring melted wax into the container through 
holes and then allowed to solidify. The coordinates of each sphere were then 
individually and painstakingly measured (probably by one of his postgraduate 
students). This seemingly ad hoc procedure yielded a maximum hard— sphere 
non— crystalline packing of 0.6366. Amazingly, more recent attemts to reproduce 
this structure using modern computational techniques have at best produced 
indistinguishable results. 

In 1968, Verlet [ref. 40] published the results of his computer experiments 
on classical fluids. His was one of the earliest attempts at modelling the liquid 
structure through a computer simulation. His model used a Lennard— Jones 
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potential as the means of particle interaction from which the equations of motion 
were solved. His study offered evidence that correlation functions at high density is 
due to the geometrical effects of a strong repulsion in the potential. He showed that 
the same behavior can be obtained by the approximate solution of the hard sphere 
problem and that the diameter of the hard spheres is the only parameter of the 
theory (in agreement with the simplistic Bernal model). 

In 1982, Miranda and Torra [ref. 41] obtained good results for the liquid 
structure and for the self-diffusion constant with a computer simulation using 
functions derived from a local pseudopotential and various other dielectric functions. 

3. Computer Simulations Summary of Liquid Sputtering 

In 1985, Lo et al [ref. 42] generated two liquid targets consisting of 603 
copper atoms and then subjected them to an argon ion beam. The two targets 
differed only in the imposed boundary conditions. One target used a box boundary 
condition requiring particles to experience pure reflection at the boundaries and the 
other used a semi— periodic boundary condition requiring position and momentum 
periodicity in the two dimensions defined by the surface. They found that the 
semi— periodic boundary condition results were in better agreement with 
experimental sputtering results and concluded that such a target better represented 
the free surface of a real liquid. The target was warmed to liquid temperatures by 
assigning velocities to each atom with a random number generator. The target was 
then allowed to equilibrate for a few picoseconds. This energy imput scheme is 
particularly important to this research effort because it is used in generating one of 
our targets. 
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In 1987, Lo et al [ref. 43] again conducted a computer simulation study of 
collision cascades in liquid indium. Their results suggested that the detailed 
structure of the target surface layer is very important in the sputtering proccess. 

The most recent simulation of liquid sputtering was done by Morgan 
[ref. 44]. He conducted a study of self— sputtering using a stratified liquid metal 
surface as a model. He observed an enhanced low energy yield which fell below 
those of other models for higher energies. These results appear to be in general 
agreement with various published measurements of liquid sputtering yields 
[refs. 34-36]. 
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II. OBJECTIVES 



This research and computer simulation effort has three specific objectives. The 
first objective is to develop a reasonable model of copper in the liquid phase by 
using a heavily modified version of QDYN (Quick Dynamics), a molecular dynamics 
computer program developed by Harrison [ref. 45]. The general intent, is to warm 
an fee (copper) crystal by giving each atom in the crystal sufficient energy to reach 
liquid temperatures. The energy given to each atom is determined stochastically. 
The warmed ensemble is then allowed to interact through pair-wise forces until 
thermal equilibrium is reached at the desired liquid temperature. The resulting 
liquid targets will be compared to various known liquid— like characteristics such as 
the radial distribution function for copper in the liquid state and the expected 
Maxwell— Boltzmann (velocity distribution) function. 

Eight different liquid targets will be generated using two different methods of 
warming (adding energy to the atoms) in combination with four different variations 
of dynamic integration. These variations are the permuted combinations of 
simulation runs with or without a reflective boundary condition and, or, with or 
without the requirement to update the nearest neighbor list in each time step. 

The second objective is to compare the dynamic behavior of the two different 
warming methods. One of these warming methods is a variation of the method used 
by Lo et al [ref. 42], where the energy is added to each atom by assigning each 
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a random velocity. The second warming method was one suggested by the late 
professor Harrison of the Naval Postgraduate School. This involves warming a 
target by randomly displacing each atom from their stable positions and then 
allowing then to interact through pair-wise forces until equilibration is reached. 
Clearly, each method starts the molecular dynamics with a different form of excess 
energy. One method starts with excess kinetic energy and the other starts with 
excess potential energy. Therefore, the ability of QDYN to equilibrate the targets 
from such varying initial energy conditions will serve as further evidence for the 
applicability of pair— potentials in molecular dynamics. 

The third objective is to use the resulting liquid targets in a computer 
sputtering experiment. The goal is to obtain sputtering yields which are in general 
agreement with published experimental data [refs. 33—36] as well as other liquid 
sputtering simulation results [refs. 42—44]. 
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III. THEORY 



A. LIQUIDS 

Matter manifests itself in three states. One of these is the liquid state and is 
the least understood of all three. Understanding the unique descriptive 
characteristics of a liquid, especially those that may affect the interatomic collision 
process, such as density, radial distribution, velocity distribution and dynamic 
behavior, is of great importance to this research effort. 

1. The Density of Liquids 

The essential difference between a solid and a liquid is that a liquid is a 
much more disordered state [ref. 47]. A disordered ensemble of atoms requires more 
space than does a typical closed— packed and well-ordered crystal structure. 
Consequently, the density of a liquid structure can be expected to have lower 
density than that of a solid. In the case of copper, the density has been 
experimentally measured for a wide range of temperatures in the liquid phase and is 
found to linearly decrease with increasing temperatures. Cahill and Kisherbaum 
[ref. 48] have obtained experimental density measurements for copper from the 
melting point (1356 K) to 2500 K. Their results indicate that a plot of the density 
of copper versus temperature in the liquid phase can be accurately described by a 
straight line of negative slope given by, 

p x = 9.077 - (8.006xl0' 4 T t ) 

where T, is the liquid temperature of copper in degrees Kelvin and p x is the density 
of copper in grams per cubic centimeter. 
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2. The Radial Distribution of Liquids 

The liquid state imposes spatial restrictions on the possible ways that 
atoms can arrange to form amorphous structures. The nearest neighbors to a given 
atom cannot move too far away because of the repulsion force exerted by the other 
particles in the system. Due to these restrictions the liquid state has, on average, a 
very distinctive radial distribution as described by Azaroff [ref. 49]. This short 
range order of liquid makes it possible to describe atomic structure in terms of the 
average density of atoms per unit volume p Q and the actual density p( r) a distance r 
away from an atom placed at the origin. 

In monatomic liquids such as copper, the radial density is spherically 
symmetric and can be described as the number of atoms per unit volume contained 
in a spherical shell of radius r and of thickness dr around one atom designated as the 
origin. The volume of such a shell is the difference between the volumes of tw r o 
spheres whose radii are r and r+dr. The density of atoms in such a shell is given by 
47rr 2 p(r)dr. 

The radial distribution function for the atoms is defined as, 

G(r) = 47rr 2 p(r). 

The radial distribution, however, cannot be experimentally measured. 
Instead, the structure factor i(s), defined below, is obtained through x-ray or 
neutron diffraction experiments from which the radial distribution is calculated. 

The radial distribution function G(r) is related to the structure factor i(s) 
by the following Fourier transformation [ref. 50]: 

9 r /' 00 

G(r) = 47rr 2 p(r) = 4xr 2 p Q + f s i(s) sin(rs) ds 

o 

where p Q is the average density of the ensemble. 
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The problem with this transformation is that the values for i(s) must be 
known through oo . To avoid this problem, the method of Zei and Steffan [ref. 51] 
can be used to extend the range of measured i(s) without manipulating the 
experimental values. This is the method used by Eder et al [ref. 46] to obtain the 
radial distribution of liquid copper at 1393K and 1S33K. This is the radial 
distribution against which the results of this simulation will be compared. 

In computer simulations, the radial distribution is very easily obtained by 
tallying the actual separations of each atom from every other atom in the ensemble. 
This is the method used in our simulation. The distribution is normalized by 
dividing each bin total by the area of a shell whose radius is equal to the minimum 
distance of the bin. The reason for doing this is to better compare the distribution 
of a very small ensemble with that of a real system whose size is infinitely larger by 
comparison. 

3. The Velocity Distribution of Liquids 

The general collision behavior of liquid particles as they interact through 
pair-wise forces is not very different from that of a gas. Therefore, we can expect a 
velocity distribution of a liquid to be similar to the velocity distribution of a gas. 

The distribution of molecular speeds in a large sample of gas varies over a 
wide range of magnitude and has a characteristic distribution which depends on 
temperature. The first expression for the velocity distribution of a gas was derived 
out by Clerk Maxwell. According to Maxwell, a sample of gas containing N 
molecules has a distribution of velocities given by [ref. 52]: 



N(v) = 47rN t 



m/27rkT 



v 2 e -(mv 2 /2kT) 



where N(v)dv is the number of molecules having speeds between v and v+dv, T is 
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the absolute temperature, k is the Boltzmann's constant, N t is the total number of 
molecules and m is the mass of the molecule. 

4. Equilibrium and the Equipartition of Energy 

The equipartition of energy theorem states that when the number of 
particles is large and Newtonian mechanics hold, all the energy terms have the same 
average value, and that the average value depends only on the temperature [ref. 52]. 
This means that the available energy distributes itself in equal shares to each of the 
independent ways in which molecules can absorb energy. 

The energy given to the atoms by the warming scheme in the simulation 
will distribute itself into half kinetic energy and half potential energy. The energy 
going into kinetic will further distribute itself equally among the x, y and z 
components of kinetic energy. The energy is distributed through the many collisions 
that the atoms experience as the ensemble comes to equilibrium. The target is at 
equilibrium when the total added energy is completely partitioned. 

A plot of the total kinetic energy versus time should resemble a damped 
sinusoid, oscillating about a value equal to half of the total added energy. A plot of 
the x,y and z components of kinetic versus time should show the components 
becoming equal as time progresses and their value approaches one sixth of the total 
added energy. 

B. SPUTTERING 

1. General Concepts 

Sputtering is the erosion of a target surface through the removal of atoms 
by the action of incident energetic particles. Sputtering occurs in nature when a hot 
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plasma, such as that found in a lightening bolt, comes in contact with a solid 
surface. Sputtering can be produced in the laboratory by subjecting a target surface 
to a plasma or a particle ion beam. 

The most obvious measure of erosion effects is through the sputtering 
yield, Y, defined as the average number of atoms removed per incident particle. 
The yield is a function of target structure, target atom mass, incident beam 
alignment (relative to target structure), energy of incident particles and the mass of 
the incident particle. Sputtered particles can leave the target's surface with a broad 
distribution of energy, charge state and exit angles. In order to obtain reproducible 
experiments, the following conditions must be met [ref. 53]: 

• The target surface must be clean. 

• The gas pressure must be low enough such that the mean free path of ions 
and sputtered atoms is large. 

• The ion current density must be high and the background pressure low so 
that formation of surface layers is prevented during the experiment. 

• The ions must strike the target at a known angle. 

• The energy spread of the incident beam must be small. 

• The ionizing conditions in the ion source should minimize the production of 
multiple charged species; the atoms must be uniformly charged and 
separated. 

• The lattice orientation of monocrystalline targets or texture 
of polycrystal must be known. 

Sputtering theory has four major variants; collisional, thermal, electronic, 
and exfoliational. In this research effort, we are primarily interested in the 
collisional behavior of sputtering. 
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2. Collisional Description 



The collisional description of sputtering views the bombarding ions as 
energetic projectiles which deposit kinetic energy and momentum upon striking a 
target atom. This action is analogous to an atomic size game of billiards. 

There are two main collisional processes; prompt collisional and slow 
collisional. Their name arises from their respective collision— event time duration. 
Once again we limit the scope of the research by concentrating on the prompt 
collisional process. 

The prompt collisional process has a duration time of approximately 500 
femptoseconds following impact. The general process occurs in the following 
manner; an incoming particle collides with a target atom and transfers some of its 
kinetic energy to it. If the energy received by the target atom exceeds the binding 
energy of its lattice site, a primary knock-on atom is created. This knock-on atom 
can travel through the target's lattice colliding with other atoms and transferring 
some of its kinetic energy and momentum. A near surface atom is sputtered if its 
kinetic energy has a component normal to the surface which is larger than the 
surface potential energy barrier. 

There are four basic collisional sputtering descriptions; rebound, recoil, 
reflection, and direct or reflected recoil, these are summarized in turn in the 
following sub— sections. 

a. Rebound Sputtering 

In the rebound sputtering sequence, an incident ion strikes an atom 
in the first atom layer. The atom initially receives a component of kinetic energy 
normal to the target surface and into the target. The atom then collides and recoils 
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off an atom in the second layer reversing the surface normal component of kinetic 
energy. The atom then escapes the surface through the first atom layer as 
illustrated by Figure 1 of Appendix A. 

This mechanism applies only to adatoms which are much lighter 
than atoms of the substrate and assumes that atoms recoil without energy loss. 

b. Recoil Sputtering 

In the recoil sputtering process, the first atom struck by the 
bombarding particle is not the ejected atom. The bombarding particle imparts 
kinetic energy and momentum to an atom in the first atom layer in a collision 
process similar to the one described in rebound sputtering. The atom then 
undergoes a sequence of collisions with other target atoms as it plows inwards. The 
atom's initial inward component of kinetic energy is eventually turned back towards 
the surface through multiple recoil collisions. The atom finally strikes another atom 
from below causing its ejection from the surface. Often, only two collisions are 
required as illustrated by Figure 2 of Appendix A. 

c. Reflection Sputtering 

In reflection sputtering, the incoming ion is backscattered by a 
target atom below layer one. The reflected ion strikes an atom in layer one from 
below causing its ejection (see Figure 3 in Appendix A). This process is extremely 
rare due to the space and angle restrictions imposed by the lattice geometry. 

d. Sputtering by Direct or Deflected Recoil 

In direct or deflected recoil sputtering, an atom in any near-surface 
layer is struck by the incoming particle at grazing incidence. The struck atom is 
deflected by neighboring atoms, causing its ejection from the target surface. This 
process is illustrated by Figure 4 of Appendix A. 
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The four collisional sputtering mechanisms described above are 
simplified and highly idealized. What most often occurs in a sputtering sequence, 
however, is the formation of collisional cascades. When the first atom in a sequence 
is struck by the bombarding particle, it will undergo a series of collisions throughout 
the crystal structure until the atom exhausts its kinetic energy or it is ejected from 
the target structure. These multiple collision sequences allow many paths through 
which the initial ion energy can reach the surface of the target and cause sputtering 
or reach deep target atom layers and dissipate its energy as heat. 

As a collision cascade develops, it may cause the sputtering of atoms 
through any of the mechanisms described above or a combination thereof. However, 
it is also possible for a cascade to produce no ejections at all. Such a cascade is 
illustrated by Figure 5 of Appendix A. 

C. THE COMPUTER MODELS 
1. Experimental Approach 
a. General Description 

The general approach to this study is to first obtain several 
reasonable liquid targets and then use the the best resulting target in a sputtering 
run. The first step is to generate a fcc(lll) crystal lattice. This is done with the 
subroutine Fill [ref. 54]. Next, the atom ensemble is energized using a stochastic 
scheme. In this simulation, there are two distinct methods used to randomly assign 
energy to individual atoms. These will be individually discussed in later sections. 

After the excess energy is given to the target, the dynamic part of 
the liquid simulation starts, and the atoms are allowed to interact through 
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pair— wise forces until they reach thermal equilibrium. This part of the simulation is 
done with a modified version of QDYN [ref. 45]. In sputtering simulations using 
QDYN, atoms are only assumed to move after collision with a moving atom, in 
generating the liquid all the atoms are designated as moving atoms from time— step 
one. 

b. Assumptions 

The most significant assumption is that the motion of particles in a 
liquid state can be approximated by pair— potential derived forces and that the 
particles obey Newtonian Mechanics. Another assumption (for sputtering) is that 
the atoms in a liquid target move so slow, compared to the ion and collision cascade 
atoms, that they can be considered static. The average velocity of liquid atoms at 
1700K is about 670 m/sec. The average collision cascade atom is at least 10 times 
faster and the ion is almost 3000 times faster. A one KeV Argon ion moves at 3030 
km/sec. So during a typical cascade lasting 300 fs the average liquid atom moves 
2xl0" 2 m, approximately one half of a lattice unit. 

2. Liquid Target Generation 

The basic target consists of 1445 atoms originally placed at the lattice 
positions of an fee crystal with a 111 plane orientation. The atom positions are 
defined in terms of lattice units. The lattice unit for copper is consistently used 
throughout the simulation as the basic unit of length. This unit, as well as other 
copper constants used in the simulation, can be found in Table 1 of Appendix B. 
The dimensions of the crystal are 19 x 8 x 17 lattice units. 
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The target is aligned in such a way that an incoming ion would strike the 
target in the 111 plane and travel in the positive y direction as illustrated by 
Figure 6 of Appendix A. 

The atoms are warmed using two techniques as explained in the following 
sub— sections. 

a. Warming by Impulse 

This method of warming atoms is a variant of the method used 
by Lo et al [ref. 42]. The general idea is to give each atom an initial impulse whose 
direction and magnitude is determined by a random number generator. 

The warming is accomplished in two steps and each uses a 
different random number generator. The magnitude of velocity for each atom is 
obtained from the product of a multiplicative constant (the speed of interest) and a 
set of random numbers having a Gaussian distribution about unity with a standard 
deviation of one. The random numbers were scaled to go from zero to two. This 
was done to avoid large velocities. The multiplicative velocity constant is chosen to 
be the velocity corresponding to twice the desired kinetic energy. This is done to 
allow for the equipartition of half of the added energy into potential energy 
elements. In our case, the desired temperature is 1462 K. Therefore, the velocity 
corresponding to twice the temperature is 1067.7 meters per second. This is the 
multiplicative constant used in the generation of random speeds. The individual 
atom velocities are the products of this speed with a series of random numbers. 

The direction of each atom must be uniformly distributed in 3— D 
space. The velocity directions are obtained from the sines and cosines of a spherical 
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coordinate system. The sines and cosines are obtained from a random number 
generator which produces uniform random numbers between zero and one. 

The actual warming is done by a subroutine called "warm", excerpt 
of which is shown below, 



SUBROUTINE WARM 
PARAMETER(LTMX=1854) 

DIMENSION RAN(LTMX),VRAN(LTMX*5) 



DO 300 1=1, LT 
R1=VRAN(I) 

R2=VRAN(I+(2*LT)) 

R3=VRAN 1+ 3*LT) 

R4=VRAN(I+(4*LT)j 

R5=VRAN(I+(5*LT)) 

VELMAG=VMAG*RAN(I) 

SIGN1=+1.0D0 

IF(R5.LT.0.5D0)SIGN1=— 1.0D0 
SIGN2=+1.0D0 

IF(R5.LT.0.5D0)SIGN2=— 1.0D0 
COSPHI=SIGN2*Rl 

SINPHI=DSQRT(1.0D0— COSPHI*COSPHI) 
SINTHE=SIGN1*(2.0D0*R3*R4)/(R3*R3+R4*R4) 
COSTHE=(R3*R3— R4*R4)/(R3*R3+R4*R4) 
VX(I)=VELMAG*(SINPHI*COSTHE) 
VY(I)=VELMAG*(SINPHI*SINTHE) 
VZ(I)=VELMAG*COSPHI 
300 CONTINUE 



CONTINUE 

RETURN 

This do-loop goes from 1 to LT, where LT is the total number of 
atoms. VELMAG is the velocity magnitude of each atom I, obtained from the 
product of a velocity multiplicative constant, VMAG, and a random number 
RAN(I). The vector RAN(I) is a set of random numbers with a gaussian 
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distribution about unity and with a standard deviation of one. These numbers are 
not allowed to be greater than two or less than zero to avoid large velocities. 
Rl— R5 are uniformly distributed sets of random numbers obtained from the vector 
VRAN(I). R2 and R5 are used to define the sign(+/— ) for SINTHE and COSPHI. 
COSPHI, SINPHI, SINTHE, and COSTHE are the sines and cosines of the usual 
spherical coordinate directional angles and a. VX(I), VY(I), and VZ(I) are the 
resulting components of velocity for atom I. 

The liquid simulation program using this type of warming technique 

is called QLV. 

b. Warming by Displacement 

This warming scheme gives added potential energy to the ensemble 
by individually displacing each atom away from their positions of minimum 
potential. The displacement is assigned to each atom by adding multiples of the 
thermal amplitude to the x,y and z components of position. The individual 
component displacement is the product of the thermal amplitude and a random 
number. The random number generator used produces numbers which are Gaussian 
distributed about zero and with a standard deviation of one. 

The thermal amplitude, T amp , is obtained from the relation, 



where u - is the mean square of total displacements of an atom from the average 
position. This is a value which can be experimentally obtained. 

The existing experimental measurements of u >l are for temperatures 
well below the liquid phase. For this reason, the warming must be done in cycles. 
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Each atom is displaced several times until a suitable liquid temperature is achieved. 
The first warming cycle consists of displacing all atoms from their equilibrium 
positions by adding to each component of position a random multiple of a known 
thermal amplitude. The following warming cycles further displaces each atom from 
their last position using the same principle. The warming behavior of such 
sequential cycling is that a maximum average temperature is reached after several 
cycles. Further displacements cause the total energy to oscillate about some 
average. 

The thermal amplitude corresponding to 300K was chosen because it 
was an actual data point in the experimental results of Singh and Sharma [ref. 55]. 
The corresponding thermal amplitude used in the simulation is 0.0782 lattice units. 
The liquid simulation program using this type of warming technique is called QLP. 
c. Achieving Thermal Equilibrium 

Thermal equilibrium is achieved by allowing the atoms to interact 
with each other through pair-wise forces. The target is considered to be in thermal 
equilibrium when the plot of kinetic energy versus time settles at half of the total 
added energy and the x,y and z components settle at equal values. Since the added 
energy in the impulse-warmed program starts with all kinetic energy and the 
displacement-warmed program starts with all potential energy, the plots of kinetic 
energy versus time for QLP and QLV resemble mirror images. These plots are 
similar to damped sinusoids which approach the value of half of the total added 
energy as the target comes to equilibrium. 
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d. Boundary Conditions 

Simulations were run with and without boundary conditions. The 
boundary conditions requires the atoms to experience pure reflection at the 
boundaries. The target boundaries for reflection are the uniformly expanded solid 
crystal dimensions corresponding to a liquid density. The density of copper is 
obtained from the experimental measurements of Cahill and Kirshenbaum [ref. 48]. 

The simulation programs were also run without boundary conditions 
to investigate the extent of the overall system expansion due to the interatomic 
interactions. 

3. Sputtering Program 

The computer program used for sputtering is SDYN88A, and it is a 
variant of QDYN developed by Harrison [ref. 45]. The program uses 
multiple-interaction (MI) logic in a time-step approach. 

The following actions occur in each time— step; 

• Summation of pair-wise forces for each atom. 

• Calculation of new velocities and positions. 

• Movement of atoms to their new positions. 

• Test energy conservation. 

The forces affecting individual atoms are not computed until a collision 
occurs. This initial collision turns the atoms on (designates them as moving atoms) 
and puts them on a moving— atom list. Once in this list, the atom's interactions are 
computed each time step. 

The moving atoms maintain a list of the atoms with which they can 
interact. This list is called the nearest neighbor list. The process of updating this 



28 



list requires considerable computing time. For this reason, it is updated every five 
time steps. 

The forces are computed by solving Newton's equations of motion. A 
predictor— corrector integration scheme adjusts the time increment of each time step 
based on the fastest atom. The integration technique is unique in that position and 
velocity are calculated for the same instant of time. 

The program was originally designed to study the effects of sputtering in 
solid targets. Consequently, the original logic did not calculate the new velocities 
and positions of every atom each time— step. Instead, only the moving atoms were 
updated. In our version, however, all of the target atoms have velocity from the 
start of dynamic integration. The initial atom positions and velocities are read 
from an input file created by the liquid simulation programs QLP and QLV. 

In order to study the validity of using a static target in liquid sputtering 
simulations, we have allowed for the simulation to proceed with either the atoms 
turned "on" or "off". When the atoms are turned on, the simulation assumes every 
atom in the target is moving. When the atoms are turned off, the individual atom 
velocities are zeroed until they undergo a collision with a fast atom or ion. Such 
static targets are the equivalent of an amorphous solid with liquid densities. The 
computing time for a run with the atoms turned on is about five times greater. The 
justification for using a static target is that the velocity of a IKeV ion is much 
faster than that of an atom in the liquid state (about 3000 times faster for 1700 K). 
The average velocity of a Cu atom in the liquid state is about 800 m/sec. The 
collision cascade atoms are about 10 times faster than liquid atoms. 
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The ion impact points are determined by a grid of 300 irregularly spaced 
points. The grid is one lattice unit squared. The grid is overlaid near the center of 
the target. For our target the grid reference point was chosen as x=6.25 LU and 
y=11.5 LU. The impact point coordinates are read in from the input file and added 
to the reference position. 

D. POTENTIAL FUNCTIONS 
1. Background 

If we take two atoms originally separated by a large distance and slowly 
bring them together, the atoms will eventually start to experience an attractive 
force. This attraction can be simulated with an attractive potential function that 
has the characteristics of a well. As the particles get closer, the potential decreases 
until the point where the force changes from attractive to repulsive. This point 
corresponds to the equilibrium separation. As the pair gets closer together, the 
repulsive force quickly increases. This is known as the hard collision range. The 
potential in this range has the characteristics of a wall. 

In the simulation, the molecular dynamic interactions are caused by 
two— body forces. The basic assumption is that the forces and potential energies 
depend solely on the physical properties of the interacting particles and their 
internuclear separation. The individual pair-wise interactions that each atom 
experiences due to its neighbors is summed to obtain the resultant force for that 
particle. 

The range of the attractive force for Cu— Cu is assumed to be 2.4 lattice 
units. This range corresponds to a point between the second and third nearest 
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neighbors. This range cut-off has been investigated and found to be reasonable for 
fee targets [ref. 56]. The differences in potential beyond the second nearest 
neighbors were found to be minor. 

2. General Forms 

The atom— atom potentials are composite functions consisting of a 
repulsive section, the wall , joined smoothly to an attractive section, the well. The 
joining of the two functions is accomplished by a cubic spline function. The 
atom— ion potential is a strictly repulsive wall, which is much steeper than that of 
the atom— atom. 






where A is the atom's hardness and B is the atom's size. 



The atom— ion repulsive wall is a Moliere potential function given by, 
V( r ) = [(ZiZ 2 e 2 /a)/(r/a)]^(r/a) 



where 



• ip{ r/a) = [0.35e~^' 3r / a + 0 55e L2r / a -f OTe -6 ’ 01- /^ 




1 - 2/3 



• a,= Bohr radius 

0 



Zi = atomic number of atom 



Z 2 = atomic number of ion. 
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The general form used for the well is a Morse potential of the form, 
V(r) = DJ e -<2a ( r-r e) > _ 2e" <a ( r - r e) > ] 

where 



• D e = well depth 

• r e = equilibrium separation 

• a = well width control. 



1. The Functions Used in the Simulation 

Idle atom— atom potential is a repulsive Born— Mayer joined to an 
attractive Morse potential by a cubic spline. 

The composite pair potential function is given by, 



Vjj — C Q +Cir -+■ C 2 r~ + C 3 r^ 

V.. = D e [ e - <2Q ( r " r e)> _2e- < ^ r " r e) > ] 



r < R a 

R a < r < R b 

R b ^ r < R c 
R c <r 



where Vjj is the composite potential for the i th and j tl ‘ atom separated by a 
distance r. 

The cubic spline joins the Born— Mayer and Morse at points R a and R b . 
The attractive Morse function is truncated at R c , between the second and third 
nearest neighbors. 



The atom— ion potential is obtained with the purely repulsive Moliere 



function, 

Vi = [(Z,Z, e 2 /a) / (r/a)] t>(r/a) r < R, 

V, = 0 r > R a 

where Vj is the potential of the i th atom and the ion. 
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The potential function parameters used in the simulation are found in 
Table 2 of Appendix B. The cubic spline parameters are found in Table 3 of 
Appendix B. The actual Ar— Cu and Cu— Cu potential functions are plotted in 
Figures 7 and 8 of Appendix A. 

E. ERROR ANALYSIS 

1. The Experimental Error 

Since the system is non dissipative, energy is conserved and the energy 
calculation results may be used as a possible check of numerical error. In sputtering 
simulations, a small percent of energy error (up to 3%) is considered acceptable. 
The justification is that studies have shown that a small change in the initial ion 
energy does not significantly affect the sputtering yield. In the simulation of a 
liquid, however, the conservation of energy is an important factor which should be 
maintained. 

The principal reason for our concern with energy conservation is that an 
error in energy is related to an error in temperature by, 

energy error = 3/2Nk x temperature error 
where N is the number of particles and k is the Boltzmann constant. 

This equation says that a 1 eV uncertainty in energy corresponds to a 
maximum uncertainty in temperature of 5.35 K. This means that a 1% error in the 
total energy of a liquid target, say 2461 eV, equilibrated to a temperature of 1575 K 
would correspond to an uncertainty in temperature of ± 65.95 K. Such an 
uncertainty in temperature would be critical if the temperature of the liquid is close 
to either the melting or boiling points. The largest energy loss in the simulation 
was 6.45 eV with a corresponding temperature uncertainty of 34.51 K. 
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2. Establishing the Fit of Curves 

According to Maxwell's equipartition of energy, the energy added to a 
number of molecules will be equally partitioned between the potential energy and 
the components of kinetic energy as the ensemble comes to thermal equilibrium. 
The variation of total kinetic energy with time is, therefore, a good indicator of 
thermal equilibrium. In our simulation, when the total kinetic energy approaches a 
value equal to half of the total added energy and when the components of kinetic 
energy become more or less equal, the group of atoms will be considered to be in 
thermal equilibrium. How close these relationships come to their theoretically 
expected values will be the principal measure of target equilibration. 

The third type of relationship generated from the results of the liquid 
simulation programs is the radial distribution of the atoms. This has been 
determined experimentally [ref. 46]. The calculated radial distribution function will 
be compared to the experimental radial distribution function by the general shape 
and by the coincidence of their peaks and valleys. 

The velocity distribution can be compared to the theoretical expected 
value using statistical analysis. The expected distribution of the atom velocities is 
the Maxwell-Boltzmann distribution. \ 2 tests will be used to determine how well 
the calculated velocity distributions fit the expected theoretical values of the 
Maxwell-Boltzmann distribution function. The following equation is used to obtain 
the x 2 values, 

(O k - E k ) 2 

El 



A 2 = 1/d X 

k = 1 
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where 



• d = number of degrees of freedom 

• k = 1 • • *n, individual number of bins 

• Ok = observed values 

• Ek = expected values. 

Generally, if a 2 » U the observed results do not fit the assumed 
distribution very well. The fit is generally considered acceptable for a a 2 < 1- 

A more qualitative measurement of the fit is made by finding the 
percentage probability Pd(A /2 > A 2 ) °f obtaining a value of a 2 > A'q- Where a 2 is the 
actual value of \ 2 obtained from experimental measurements. These probabilities 
are found in Table D, page 251 of reference 58. If the probability Pd(x 2 > A 2 ) < 5% 
then the fit is considered to be acceptable. This is the criteria applied to our 
results. 

In our simulation, the number of bins used for the velocity distribution is 
30 and the degrees of freedom is d = 29. Using Table D of reference 58, the 
maximum a / 2 corresponding to a Pd(A /2 > A~) = 5% for 29 degrees of freedom is 
1.493. Therefore, if the a 2 of any of our velocity distributions is under 1.493, then 
the fit will be considered to be consistent with the expected values. 
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IV. RESULTS AND DISCUSSION 



Eight different targets were generated using various boundary conditions, 
warming methods and nearest neighbor list update schemes. The target with the 
best liquid— like characteristics was chosen and subjected to a 295 trajectory 
sputtering run. 

In order to facilitate reference to the various simulation schemes and their 
results the following names are defined, 



• QLP, Liquid simulation using the atom displacement warming scheme and 
unrestrained boundaries. 

• QLPBC, Liquid simulation using the atom displacement warming scheme 
and reflective boundary conditions. 

• QLP— REVl, Liquid simulation using the atom displacement warming 
scheme, unrestrained boundaries and a new neighbor— list update every 
time step. 

• QLPBC— RENT, Liquid simulation using the atom displacement warming 
scheme, reflective boundary conditions and a new neighbor-list update 
every time step. 

• QLY, Liquid simulation using the impulse warming scheme and unrestrained 
boundaries. 

• QLYBC, Liquid simulation using the impulse warming scheme 
and reflective boundary conditions. 

• QLY— REV1, Liquid simulation using the impulse warming scheme, 
unrestrained boundaries and a new neighbor— list update every time step. 

• QLVBC— RE VI, Liquid simulation using the impulse warming scheme, 
reflective boundary conditions and a new neighbor— list update 

every time step. 
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BC means that reflective boundary conditions are applied and REV1 stands for 
the revision which includes the update of the nearest neighbor list every time step. 

A. ENERGY CONSERVATION 

The calculated total energy losses for all eight runs range from 0.001% to 0.25%. 
The equivalent energy uncertainty for this range is 0.03 eV to 6.45 eV. Although all 
eight results have acceptable energy losses, there is a clear division in the energy 
conservation efficiency betweeen the runs that update the nearest neighbor list every 
time step and those that do not. The simulation runs which update the nearest 
neighbor list each time step give the poorest energy conservation. This is to be 
expected because as the simulation proceeds, atoms find new neighbors with which 
they can interact and this requires additional calculations that add to the total 
error. 

The best energy conservation was obtained with QLPBC. A summary of the 
experimental error for all eight targets is found in Table 4 of Appendix B. 

B. DENSITY RESULTS 

The density was calculated for each of the unrestrained targets (without 
reflective boundary conditions) after equilibration was reached. These were found 
to be slightly lower than those corresponding to liquid copper at the equilibration 
temperature. 

The most accurate density was obtained with QLV— REV1. This simulation 
yielded a liquid density of 7.664 gm/cm 3 , within 1.8% of the expected value for a 
liquid temperature of 1690 I\. Table 5 of Appendix B shows the resulting densities 
for all runs with unrestricted boundaries and their expected theoretical values. 
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C. TARGET EQUILIBRATION AND FINAL TEMPERATURES 

The plots of total kinetic energy versus elapsed time (Figs. 9—16, Appendix A) 
and those of the components of kinetic energy versus elapsed time 
(Figs. 17—24, Appendix A) show that all targets reach equilibrium after 
approximately 400 femptoseconds. No significant change in the kinetic energy or in 
the components of kinetic energy is observed after this time. This equilibration 
time seems to be consistent with the simulation results of Lo et al [ref. 42], The 
simulations were computed for at least twice the equilibration time. 

These kinetic energy plots also show that equipartition of energy holds very well 
for all eight target generation schemes. In particular, QLPBC produced a final 
kinetic energy of 273 ± 17 eV, less than 1 eV difference from the theoretical expected 
value of 273.071. The error for the total kinetic energies is taken from the 
maximum amplitude of the kinetic energy oscillations about an average value 
following the equilibration point. Table 6 of Appendix B summarizes the final 
kinetic energy and temperature results for all eight targets. 

QLV-REY1 and QLVBC-REVl produce the smallest kinetic energy 
oscillations at equilibrium, ±4 eV. This corresponds to a temperature uncertainty of 
±32 K. The final temperatures for these two targets are 1690 ± 32 I< and 1700 ± 32 
respectively. The reason for such oscillations in kinetic energy after equilibration is 
reached is that the size of the ensemble is small compared to a real liquid system. 
Real systems have many billions of particles and the total kinetic energy 
fluctuations can average out. 
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D. RADIAL DISTRIBUTIONS 

The radial distribution for an unwarmed solid copper crystal is given by 
Figure 25 of Appendix A. A comparison of this figure with the various radial 
distributions produced by the liquid simulations (Figs. 26—33 of Appendix A) 
clearly contrasts the structure differences between a solid and a liquid. These liquid 
radial distributions have the distinctive wide and smooth, peaks and valleys 
characteristic of a liquid. The crystal structure is completely lost in all samples. 

The peaks and valleys of all the resulting radial distributions are compared 
against each other for consistancy, and against neutron diffraction experimental 
results [ref. 46] for validation, in Table 7 of Appendix B. All the numbers for this 
table are taken directly from the radial distributions plots, Figures 26—33 of 
Appendix A. It is important to note that many of the expected valleys could not be 
positively identified and were therefore left blank. 

The best (most liquid— like) radial distributions are those obtained with 
QLP-REV1, QLPBC-REV1, QLV-REV1, and QLVBC-REV1. These radial 
distributions are in excellent agreement with the neutron diffraction experimental 
results [ref. 46]. The very best radial distribution of all is that of QLVBC— REV1. 

E. VELOCITY DISTRIBUTIONS 

The resulting velocity distribution and the expected Maxwellian distribution for 
each target are plotted together in Figures 34^41 of Appendix A. These figures 
show that the velocity distribution of all eight liquid targets generally match the 
Maxwellian distribution. The reduced x 2 tests and their probabilities PdU' 2 ^X' 2 ) are 
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found in Table 8 of Appendix B. This table shows that the most Maxwellian-like 
velocity distribution is obtained with QLP. Only three targets QLP, QLV and 
QLVBC— REV1 have velocity distributions with a Pd(A 2 ^X 2 ) > 5%. This is the 
standard rejection criteria. The smallest Pd(A /2 >X 2 ) was 0-5% with most falling 
around 1.5%. In general, the velocity distributions appear to be consistent with the 
Maxwellian distribution. 

F. SPUTTERING RESULTS 

The QLVBC— REVl target was chosen for the sputtering run because it has the 
best overall liquid — like characteristics. 295 trajectories were run with 1.0 KeV 
argon ions. The liquid target atom velocities were zeroed to save both 
computational costs and real time running requirements. A sputtering run with live 
atoms (all atoms initially moving) takes on average about five times longer than a 
static target. 

The sputtering program had on average a 3% energy loss per trajectory. 
Although this is considered acceptable for a sputtering simulation, it was later found 
that the program overestimated the amount of energy that left the target with every 
ejected atom, but did not affect the forces. Thus 3% is very much an upper bound. 

The final average yield for the 295 trajectories is 8.28 atoms per incident ion. 
This is 40% higher than the typical yield expected from a solid target on the 010 
face and at room temperature. This increase in yield is in qualitative agreement 
with experimental results [refs. 33—35] and with another computer simulation 
[ref. 43]. Table 9 in Appendix B shows the general consistancy of these results. 
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V. CONCLUSIONS 



Results from the various liquid target generation schemes have shown that the 
warming method used to energize the individual target atoms does not affect the 
final liquid-like characteristics of the resulting target. That is, results of liquid 
targets warmed by impulse were essentially indistinguishable in radial distribution 
and in velocity distribution from those warmed by displacement. Both warming 
methods are equally capable of generating liquid targets. 

The best results were obtained with the application of reflection boundary 
conditions in conjunction with a scheme which updated the nearest neighbor list 
every time step. QLPBC-REV1 and QLVBC-REV1 produced good liquid targets. 
The liquid target generated by QLVBC-REV1 is found to be the best of the two. 
Figure 43 shows the actual liquid structure generated by QLVBC-REV1. The 
straight lines in this figure represent the outline of the crystal before warming and 
equilibration. 

Results indicate that the motion of a typical liquid atom in a time span 
equivalent to at least twice the equilibration time (approximately 800 fs) is such 
that it seldom leaves its original atom neighborhood. Figure 42 in Appendix B 
illustrates this behavior. 

The sputtering of the best liquid target with a 1 KeV incident argon ion 
resulted in yields which were 40% higher than expected by solid targets. This is in 
qualitative agreement with the results of other experimental efforts [refs. 33—35] 
found in Table 9 of Appendix B. Figures 44 and 45 in Appendix A show the yield 
per incident particle distribution for both a crystal Cu(010) target and a Liquid Cu 
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target respectively. The large number of trajectories that have no yield at all for 
the crystal target may be attributed to ion channeling. 

The spot pattern (spatial distribution of ejected material) of the liquid target, 
Figure 46, shows no signs of an ordered structure. The spot pattern of a Cu(010) 
crystal, Figure 47, is shown to contrast the difference in spot pattern between a well 
ordered structure and an amorphous ensemble. 

In summary, the liquid target generation methods QLVBC— REV1 and 
QLPBC— REVl both produced reasonable liquid targets. The frequent update of the 
nearest neighbor list as well as the application of boundary conditions was found to 
be critical in making a reasonable liquid structure. The sputtering results are in 
general agreement with published data. 
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VI. RECOMMENDATIONS 



The primary purpose of this investigation was to devise a method that would 
produce a reasonable liquid target. This was achieved. 

The secondary objective was to conduct a preliminary study of liquid sputtering 
using the best liquid target. The average sputtering yield was obtained from a set 
of 295 trajectories using only one liquid target. Other targets should also be used in 
sputtering runs in order to obtain additional evidence for the validity of the liquid 
targets. It is possible that by using the same liquid target for every trajectory we 
may be introducing a systematic bias in the sputtering yield. It is reasonable to 
assume, for example, that our yields are closer to that of an amorphous solid than 
that of a liquid because the target has constant structure. 

A scheme should be devised by which several equivalent liquid targets are used 
in turn (cycled) with each trajectory of a full sputtering run. Another alternative 
would be to complete several sets of trajectories on sufficiently separated areas of 
the surface. 

At least one liquid target should be subjected to various ion energies to 
determine if there is an energy dependent yield profile as suggested by the results of 
reference 33. The effect on yield by varying the angle of incidence should also be 
investigated. 



43 



APPENDIX A - FIGURES 




Figure 1. The Rebound Sputtering Process 




Figure 2. The Recoil Sputtering Process 
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Figure 3. The Reflection Sputtering Process 
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Figure 6. The Target-Ion Alignment 
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Figure 7. Ar— Cu Potential Function. 



47 



SEPARATION (LU) 



Cu-Cu POTENTIAL FUNCTION 



cm 



r- id 




iT> 



~r 

r> 



rsi 



ID 



CM 



CO 

d 



d 



CM 



(Ae> A0d3N3 1Vi±N3±Od 

Figure 8. Cu— Cu Potential Function. 
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Figure 9. Kinetic Energy Versus Elapsed Time, QLP Output. 
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Figure 10. Kinetic Energy Versus Elapsed Time, QLPBC Output. 
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Figure 11. Kinetic Energy Versus Elapsed Time, QLP— REV1 Output. 
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Figure 12. Kinetic Energy Versus Elapsed Time, QLPBC— REVl Output. 
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Figure 13. Kinetic Energy Versus Elapsed Time, QLV Output. 
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Figure 14. Kinetic Energy Versus Elapsed Time, QLVBC Output. 
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Figure 15. Kinetic Energy Versus Elapsed Time, QLV— REV1 Output. 
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Figure 16. Kinetic Energy Versus Elapsed Time, QLVBC— REV1 Output. 
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Figure 17. Components of KE Versus Elapsed Time, QLP Output. 
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Figure 18. Components of KE Versus Elapsed Time, QLPBC Output. 
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Figure 19. Components of KE Versus Elapsed Time, QLP— REV1 Output. 
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Figure 20. Components of KE Versus Elapsed Time, QLPBC— REV1 Output. 
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Figure 21. Components of KE Versus Elapsed Time, QLV Output. 
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Figure 22. Components of KE Versus Elapsed Time, QLVBC Output. 
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Figure 23. Components of KE Versus Elapsed Time, QLV-REV1 Output. 
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Figure 24. Components of KE Versus Elapsed Time, QLVBC— REV1 Output. 
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Figure 25. Radial Distribution, Unwarmed Crystal. 
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Figure 26. Radial Distribution, QLP Output. 
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Figure 27. Radial Distribution, QLPBC Output. 
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Figure 28. Radial Distribution, QLP— REVl Output. 
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Figure 29. Radial Distribution, QLPBC— REV1 Output. 
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Figure 30. Radial Distribution, QLV Output. 
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Figure 31. Radial Distribution, QLVBC Output. 
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Figure 32. Radial Distribution, QLV— REV1 Output. 
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Figure 33. Radial Distribution, QLVBC— REV1 Output. 
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Figure 34. Velocity Distribution, QLP Output. 
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Figure 35. Velocity Distribution, QLPBC Output. 
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Figure 36. Velocity Distribution, QLP— REV1 Output. 
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Figure 37. Velocity Distribution, QLPBC— REV1 Output. 
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Figure 38. Velocity Distribution, QLV Output. 
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Figure 39. Velocity Distribution, QLVBC Output. 
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Figure 40. Velocity Distribution, QLV— REV1 Output. 
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Figure 41. Velocity Distribution, QLVBC— REV1 Output. 
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Figure 42. Motion of the Center— Most Atom Under the Influence of 
its Nearest Neighbors. Output from QLVBC— REV1. 
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Figure 43. Liquid Copper (Warmed and Equilibrated), QLVBC— REV1 
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Figure 44. Ejected Atoms per Single Ion for Solid Target, Cu(010) Crystal 







Figure 45. Ejected Atoms per Single Ion for Liquid Target QLVBC-REVl 
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Figure 46. Spot Pattern for Liquid Target QLVBC— RE\T. 




Figure 47. Spot Pattern for Solid Target, Cu(010) Crystal. 
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APPENDIX B - TABLES 



TABLE 1 


PHYSICAL DATA FOR COPPER 


! Atomic Number (Z) 


29 


i Atomic Weight 


63.540 amu 


Crystal Type 


FCC 


Lattice Constant 


3.6150 angstroms 


Re 


1.S075 angstroms 


Melting Point 


1083.4 ± 0.2 oc 


Boiling Point 


2567 oc 


Density (solid) 


8.96 gm/cm 3 



Note: Data for Table 1 is derived from reference 59. 
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TABLE 2 


POTENTIAL FUNCTION PARAMETERS 


Parameter 


Cu— Cu 


Cu-Ar 


A (KeV) 


22.564 




B (A' 1 ) 


-5. OSS 




D e (eV) 


0.431 




Re (A) 


2.628 




a (A -1 ) 


1.405 




Ra (LU) 


0.83 


1.41 


Rb (LU) 


1.10 


1.41 


R c (LU) 


2.40 


1.41 


Zi (amu) 


29 


18 


Zo (amu) 


29 


29 


I< 


0.0 


0.0920 



TABLE 3 


CUBIC SPLINE PARAMETERS 


Parameter 


Cu-Cu 


Co 


5.876x102 


Cl 


-1.594xl0 3 


C 2 


1.450x103 


C 3 


-4.423xl0 2 
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TABLE 4 


TOTAL ENERGY LOSSES 


Target 


Percent Error(%) 


Energy Uncertainty(eV) 


QLP 


0.001 


0.04 


QLPBC 


< 0.001 


0.03 


QLP-REV1 


0.251 


6.44 


QLPBC-REV1 


0.238 


6.45 


QLV 


0.01 


0.30 


QLVBC 


0.01 


0.29 


QLV-REV1 


0.23 


6.22 


QLVBC-REV1 


0.24 


4.20 



1 TABLE 5 


1 DENSITY RESULTS (TARGETS WITH UNRESTRICTED BOUNDARIES) 




Density (gm/cm 3 ) 


Target 


Expected 


Result 


QLP 


7.906 


6.646 


QLP-REV1 


7.906 


6.617 


QLV 


7.806 


6.057 


QLV-REV1 


7.806 


7.664 
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TABLE 6 


KINETIC ENERGY AND TEMPERATURE RESULTS 




KE (eV) 


Temperature (I<) 


Target 


Expected 


Result 


Expected 


Result 


QLP 


273.071 


265 ± 6 


1462.043 


1420 ± 32 


QLPBC 


i 273.071 


273 ± 6 : 


1462.043 


1460 ± 32 


QLP-REV1 


273.071 


280 ± 17 


1462.043 


1520 ± 91 


QLPBC-REV1 


273.071 


290 ± 8 


1462.043 


j 1560 ± 43 


QLV 


296.45S 


288 ± 9 


1587.258 


1540 ± 48 


QLVBC 


296.458 


299 ± 12 


1587.258 


1600 ± 64 


QLV-REV1 


296.458 


315 ± 4 


1587.258 


1690 ± 21 


QLVBC-REV1 


296.458 


280 ± 4 


1587.258 


1700 ± 21 
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TABLE 7 

RADIAL DISTRIBUTIONS, PEAKS AND VALLEYS 



Target 

QLP 

QLPBC 

QLP-REVT 

QLPBC-REV1 

QLV 

QLVBC 

QLV-REV1 

QLVBC-REVl 

NEUTRON 
DIFF. DATA 



1st 


1 st 


2nd 


2nd 


3rd 


Peak 


Valley 


Peak 


Valley 


Peak 


(LU) 


(LU) 


(LU) 


(LU) 


(LU) 



2.5±.l 

2.5±.l 

2.4±.l 

2.4±.l 

2.5±.l 

2.5±.l 

2.4±.l 

2.4±.l 

2.5±.l 



3.3±.l 



3.3±.l 

3.3±.l 

3.5±.l 



4.3±.l 

4.7±.l 

4.3±.l 

4.4±.l 

4.3±.l 

4.3±.l 

4.3±.l 

4.3±.l 

4.7±.l 



5.6±.l 



5.6±.l 



6.7±.l 

6.7±.l 

6.6±.l 

6.6±.l 

6.7±.l 

6.7±.l 

6.6±.l 

6.5±.l 

6.8±.l 



[ref. 46] 



1 



3rd 

Valley 

(LUJ 



7.7±.l 

7.7±.l 

7.9±.l 
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TABLE 8 


VELOCITY DIST. 


X 2 TESTS TO MAXWELLIAN DIST. 




Target 


X 2 


Pd(* 2 >X' 2 o) 


QLP 


0.9 


61.5 


QLPBC 


1.7 


1.6 


QLP-REV1 


1.7 


1.6 


QLPBC-REV1 


1.9 


0.5 


QLV 


1.5 


5.4 


QLVBC 


1.7 


1.6 


QLV-REV1 


1.5 


5.4 


QLVBC-REV1 


1.5 


1.5 



TABLE 9 



PERCENT INCREASE IN YIELD 

OF LIQUID SPUTTERING OVER SOLID SPUTTERING 



Source 


Type of 
Study 


Target 


Ion 


Energy 

(KeV) 


Percent Yield 
Increase (%) 


Ref. 33 


experiment 


liquid Sn 


Argon 


0.2 


40 










0.4 


-6 


Ref. 34 


experiment 


liquid Sn 


Argon 


0.2 


50 










0.375 


15 


Ref. 35 


experiment 


liquid In 


Argon 


0.107 


10 


I Ref. 43 


simulation 


liquid Cu 


Argon 


5.0 


60 


This 












Study 


simulation 


liquid Cu 


Argon 


1.0 


40 
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